
# curcomp <- "C:/Users/iosgood/Dropbox/Divisions over supply chain/IFoT_rep" 
# setwd(curcomp)

######################
## Load in the data ##
######################

# load position data
pos <- read.csv("data/IFoT_allpos_data.csv")

# load trade data
trd <- read.csv("data/IFoT_alltrd_data.csv")
trd$nonrpca0509 <- 3; rel <- trd$exp0509 - trd$nonrpimp0509
trd$nonrpca0509[rel < quantile(rel, .4)] <- 2; trd$nonrpca0509[rel < quantile(rel, .2)] <- 1; 
trd$nonrpca0509[rel > quantile(rel, .6)] <- 4; trd$nonrpca0509[rel > quantile(rel, .8)] <- 5; 

dat <- cbind(pos, trd)

dat$diff <- factor(dat$diff, levels = c("Homogeneous","Mod. differentiated","Differentiated"))

dat$posindex <- (dat$numsupfirm + dat$numoppfirm)/(dat$numsupfirm + dat$numoppfirm + dat$numsupassoc + dat$numoppassoc)

dat$lobindex <- dat$numfirmlob/(dat$numfirmlob + dat$numassoclob)
dat$lobrepindex <- dat$numrepfirmlob/(dat$numrepfirmlob + dat$numrepassoclob)
dat$lobspndindex <- dat$spndfirmlob/(dat$spndfirmlob + dat$spndassoclob)
dat$korrelus <- rnorm(nrow(dat))

dat[is.nan(dat$lobindex),c("naicsdes","naicscode")]
dat <- dat[-79,]

## load in various functions and libraries
# source("code/IFoT_repcode_Rfunctions.R")

# setwd
setwd(paste(curcomp, "/paper", sep = ""))

###########
## Table ##
###########

## Table4

cur <- dat
mod0a <- glm(I(100*posindex) ~ log(rpimp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
mod0b <- glm(I(100*posindex) ~ log(totusinp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
mod0c <- glm(I(100*posindex) ~ log(rpimp0509+1) + log(totusinp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
mod0d <- glm(I(100*posindex) ~ diff + nonrpca0509 + log(sales0509+1), data = dat)
1- pchisq(summary(mod0d)$deviance - summary(mod0c)$deviance, df = 2)

R2 <- c(myround(pR2(mod0a), 2)[6], myround(pR2(mod0b), 2)[6])
mod1a <- glm(I(100*lobindex) ~ log(rpimp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
mod1b <- glm(I(100*lobindex) ~ log(totusinp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
mod1c <- glm(I(100*lobindex) ~ log(rpimp0509+1) + log(totusinp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
mod1d <- glm(I(100*lobindex) ~ diff + nonrpca0509 + log(sales0509+1), data = dat)
1- pchisq(summary(mod1d)$deviance - summary(mod1c)$deviance, df = 2)
R2 <- c(R2, myround(pR2(mod1a), 2)[6], myround(pR2(mod1b), 2)[6])

mod2a <- glm(I(100*lobrepindex) ~ log(rpimp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
mod2b <- glm(I(100*lobrepindex) ~ log(totusinp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
# R2 <- c(R2, myround(pR2(mod2a), 2)[6], myround(pR2(mod2b), 2)[6])

mod3a <- glm(I(100*lobspndindex) ~ log(rpimp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
mod3b <- glm(I(100*lobspndindex) ~ log(totusinp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
mod3c <- glm(I(100*lobspndindex) ~ log(rpimp0509+1) + log(totusinp0509+1) + diff + nonrpca0509 + log(sales0509+1), data = dat)
mod3d <- glm(I(100*lobspndindex) ~ diff + nonrpca0509 + log(sales0509+1), data = dat)
1- pchisq(summary(mod3d)$deviance - summary(mod3c)$deviance, df = 2)
R2 <- c(R2, myround(pR2(mod3a), 2)[6], myround(pR2(mod3b), 2)[6])

mods <- list(summary2(mod0a), summary2(mod0b), summary2(mod1a), summary2(mod1b), summary2(mod3a), summary2(mod3b))
vars <- c(rownames(summary2(mod3a))[2],rownames(summary2(mod3b))[2:6]) 
tab <- matrix(data = NA, ncol = length(mods), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods)){mod <- mods[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
  stars <- rep("", nrow(mod[tabinmod,])); stars[mod[tabinmod,4] < .05] <- "^{*}"; stars[mod[tabinmod,4] < .01] <- "^{**}";
  stars[mod[tabinmod,4] < .001] <- "^{***}";  
  tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],2), stars, sep = ""), paste("(", myround(mod[tabinmod,2],2), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1,3,5,7,9,11)] <- c("Related-party imports","Imported inputs",
  "Mod. differentiated", "Differentiated", "Comp. Adv.", "Sales") # , "No. assocs.", "No. firms", "4-firm concentration", "20-firm concentration", 
N <- c(sum(summary(mod0a)$df[2:3]), sum(summary(mod0b)$df[2:3]), sum(summary(mod1a)$df[2:3]), sum(summary(mod1b)$df[2:3]), 
  sum(summary(mod3a)$df[2:3]), sum(summary(mod3b)$df[2:3]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = "")
R2 <- paste("\\multicolumn{1}{c}{", R2, "}", sep = "")
tab <- rbind(tab, c(R2), c(N)); rownames(tab)[c(nrow(tab)-1,nrow(tab))] <- c("pseudo R$^2$", "N")
addtorow <- list(); addtorow$pos <- list(12); addtorow$command <- c(paste0('\\midrule '))
ptable(cbind(rownames(tab), tab), "tables/alltradeindices.tex")

cur <- dat
model = posindex ~ log(rpimp0509+1) + diff + nonrpca0509  + log(sales0509)
csvar = "rpimp0509"; csquant = c(.25, .75); cs1 <- lmcs2(model, cur, csvar, csquant); cs1 <- cs1[[1]]
csvar = "diff"; csquant = c("Homogeneous", "Mod. differentiated"); cs3a <- lmcs2(model, cur, csvar, csquant); cs3a <- cs3a[[1]]
csvar = "diff"; csquant = c("Homogeneous", "Differentiated"); cs3b <- lmcs2(model, cur, csvar, csquant); cs3b <- cs3b[[1]]

model = posindex ~ log(totusinp0509+1) + diff + nonrpca0509  + log(sales0509) # + country2
csvar = "totusinp0509"; csquant = c(.25, .75); cs2 <- lmcs2(model, cur, csvar, csquant); cs2 <- cs2[[1]]
inf1 <- c(cs1, cs2, cs3a, cs3b)

cur <- dat
model = lobindex ~ log(rpimp0509+1) + diff + nonrpca0509  + log(sales0509) # + country2
csvar = "rpimp0509"; csquant = c(.25, .75); cs1 <- lmcs2(model, cur, csvar, csquant); cs1 <- cs1[[1]]
csvar = "diff"; csquant = c("Homogeneous", "Mod. differentiated"); cs3a <- lmcs2(model, cur, csvar, csquant); cs3a <- cs3a[[1]]
csvar = "diff"; csquant = c("Homogeneous", "Differentiated"); cs3b <- lmcs2(model, cur, csvar, csquant); cs3b <- cs3b[[1]]

model = lobindex ~ log(totusinp0509+1) + diff + nonrpca0509  + log(sales0509) # + country2
csvar = "totusinp0509"; csquant = c(.25, .75); cs2 <- lmcs2(model, cur, csvar, csquant); cs2 <- cs2[[1]]
inf1 <- c(cs1, cs2, cs3a, cs3b)

model = lobindex ~ log(rpimp0509+1) + log(totusinp0509+1) + diff + nonrpca0509  + log(sales0509) # + country2
mod <- glm(model, data = cur); vcmat <- vcov(mod); 
betas <- mvrnorm(nrow(cur), mu = coef(mod), Sigma = vcmat); 
mm1 <- apply(model.matrix(mod), 2, median); mm0 <- mm1; 
mm0[c("diffMod. differentiated")] <- 0; mm0[c("diffDifferentiated")] <- 0; mm0[c("log(rpimp0509 + 1)")] <- quantile(log(cur$rpimp0509+1), .1); mm0[c("log(totusinp0509 + 1)")] <- quantile(log(cur$totusinp0509+1), .1);
mm1[c("diffMod. differentiated")] <- 0; mm1[c("diffDifferentiated")] <- 1; mm1[c("log(rpimp0509 + 1)")] <- quantile(log(cur$rpimp0509+1), .9); mm1[c("log(totusinp0509 + 1)")] <- quantile(log(cur$totusinp0509+1), .9);
sup0 <- c(); sup1 <- c()
for(i in 1:nrow(betas)){cf0 <- mm0%*%betas[i,]; sup0[i] <- median(cf0); cf1 <- mm1%*%betas[i,]; sup1[i] <- median(cf1)}
infALL <- c(myround(median(sup1), 3), myround(median(sup0), 3), myround(quantile((sup1-sup0), .5), 3), 
  paste("[", myround(quantile((sup1-sup0), .025), 3), ", ", myround(quantile((sup1-sup0), .975), 3), "]", sep = ""))


#############
## Figures ##
#############

setwd(paste(curcomp, "/paper", sep = ""))

# basic descriptives
sum(trd$rpimp0509)/sum(trd$imp0509)
sum(trd$rpimp1014)/sum(trd$imp1014)
sum(trd$rpexp0509)/sum(trd$exp0509)
sum(trd$rpexp1014)/sum(trd$exp1014)

# manu
manu <- substr(trd$naicscode,1,1) == "3"
sum(trd$rpimp0509[manu])/sum(trd$imp0509[manu])
sum(trd$rpimp1014[manu])/sum(trd$imp1014[manu])


naics3codes <- unique(substr(trd$naicscode, 1, 3))
naics3des <- c("Plant crops", "Animal production", "Oil and gas extraction","Other mining",
  "Food", "Beverage and tobacco","Yarns and fabrics","Textile products","Apparel","Leather goods",
  "Wood products", "Paper goods", "Oil and coal prods.","Chemicals","Plastics and rubber",
  "Nonmetallics products", "Primary metals", "Fabricated metals", "Machinery",
  "Computers and electronics","Electrical and appliances","Transporation equipment",
  "Furniture", "All other manufacturing")
sum3 <- data.frame(naicscode = naics3codes, naicsdes = naics3des, 
  imp0509 = NA, rpimp0509 = NA, exp0509 = NA, rpexp0509 = NA,
  imp1014 = NA, rpimp1014 = NA, exp1014 = NA, rpexp1014 = NA)

for(i in naics3codes){
  sum3$imp0509[sum3$naicscode == i] <- sum(trd$imp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$rpimp0509[sum3$naicscode == i] <- sum(trd$rpimp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$exp0509[sum3$naicscode == i] <- sum(trd$exp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$rpexp0509[sum3$naicscode == i] <- sum(trd$rpexp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$imp0509[sum3$naicscode == i] <- sum(trd$imp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$rpimp0509[sum3$naicscode == i] <- sum(trd$rpimp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$exp0509[sum3$naicscode == i] <- sum(trd$exp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$rpexp0509[sum3$naicscode == i] <- sum(trd$rpexp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
}

# pdf("figs/rptrade.pdf", width = 4, height = 3)
tiff("figs/rptrade.tiff", width = 4, height = 3, units = 'in', res = 600)

par(mar=c(1,1.1,1.1,1.1)); # c(bottom, left, top, right
plot(-1, -1, xlim = c(-.33,1), ylim = c(0,1), axes = FALSE, xlab = "Proportion", ylab = "Industry",
  main = "                           Proportion of US Imports and Exports by Related Parties", cex.main = .5, cex.lab = .5)
text(label = c(naics3des, "All tradables"), x = -.4, y = c(c(25:2)/25, 0), pos = 4, cex = .3)
axis(1, at = c(0,.2,.4,.6,.8,1), labels = c(".0",".2",".4",".6",".8","1"), cex.axis = .3, lwd =.5, padj = -6, tck = -.015)

offset <- .005
segments(0, c(c(25:2)/25,0)+offset, c(sum3$rpimp0509/sum3$imp0509, sum(sum3$rpimp0509)/sum(sum3$imp0509)), c(c(25:2)/25,0)+offset, lwd = 1, col = "black",lend =2)
segments(0, c(c(25:2)/25,0)-offset, c(sum3$rpexp0509/sum3$exp0509, sum(sum3$rpexp0509)/sum(sum3$exp0509)), c(c(25:2)/25,0)-offset, lwd = 1, col = "gray",lend =2)

legend(x=.6,y=1, lwd = c(1,1), col = c("black","gray"), 
  legend = c("Proportion US imports rel. party","Proportion US exports rel. party"), cex = .35, pt.cex = .5)

dev.off()

sum3$imp0509[sum3$naicscode == "336"]
sum3$rpimp0509[sum3$naicscode == "336"]/sum3$imp0509[sum3$naicscode == "336"]

##
names(trd)
sum3$totusinp0509 <- NA; sum3$totusnonrpinp0509 <- NA; sum3$totusrpinp0509 <- NA; 
sum3$totusinp1014 <- NA; sum3$totusnonrpinp1014 <- NA; 

for(i in naics3codes){
  sum3$totusinp0509[sum3$naicscode == i] <- sum(trd$totusinp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$totusnonrpinp0509[sum3$naicscode == i] <- sum(trd$totusnonrpinp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$totusrpinp0509[sum3$naicscode == i] <- sum(trd$totusrpinp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$sales0509[sum3$naicscode == i] <- sum(trd$sales0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$IMPANDNONIMPtotusinp0509[sum3$naicscode == i] <- sum(trd$IMPANDNONIMPtotusinp0509[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$totusinp1014[sum3$naicscode == i] <- sum(trd$totusinp1014[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$totusnonrpinp1014[sum3$naicscode == i] <- sum(trd$totusnonrpinp1014[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$sales1014[sum3$naicscode == i] <- sum(trd$sales1014[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
  sum3$IMPANDNONIMPtotusinp1014[sum3$naicscode == i] <- sum(trd$IMPANDNONIMPtotusinp1014[substr(trd$naicscode,1,3) == i], na.rm = TRUE)
}


# pdf("figs/inttrade.pdf", width = 4, height = 3)
tiff("figs/inttrade.tiff", width = 4, height = 3, units = 'in', res = 600)

par(mar=c(1,1.1,1.1,0.5)); # c(bottom, left, top, right
plot(-1, -1, xlim = c(-.1,.4), ylim = c(0,1), axes = FALSE, xlab = "Proportion", ylab = "Industry",
  main = "                           Proportion of US Inputs Imported from Abroad", cex.main = .5, cex.lab = .5)
text(label = c(naics3des, "All tradables"), x = -.125, y = c(c(25:2)/25, 0), pos = 4, cex = .3)
axis(1, at = c(0,.1,.2,.3,.4), labels = c(".0",".1",".2",".3",".4"), cex.axis = .3, lwd =.5, padj = -6, tck = -.015)

offset <- .005
segments(0, c(c(25:2)/25,0)+offset, c(sum3$totusinp0509/sum3$sales0509, sum(sum3$totusinp0509)/sum(sum3$sales0509)), c(c(25:2)/25,0)+offset, lwd = 1, col = "black", lend = 2)
segments(0, c(c(25:2)/25,0), c(sum3$totusinp0509/sum3$IMPANDNONIMPtotusinp0509, sum(sum3$totusinp0509)/sum(sum3$IMPANDNONIMPtotusinp0509)), c(c(25:2)/25,0), lwd = 1, col = "gray", lend = 2)
segments(0, c(c(25:2)/25,0)-offset, c(sum3$totusrpinp0509/sum3$IMPANDNONIMPtotusinp0509, sum(sum3$totusrpinp0509)/sum(sum3$IMPANDNONIMPtotusinp0509)), c(c(25:2)/25,0)-offset, lwd = 1, col = "gray", lty = "dashed", lend = 2)
segments(x0 = .44, x1 = .40, y0 = (25-12)/25, y1= (25-12)/25, lwd = 1, col= "white")
text(.41, (25-12)/25, ".48", col = "gray", cex = .3)

legend(x=.225,y=1, lwd = c(1,1,1), col = c("black","gray","gray"), lty=c(1,1,2),
  legend = c("Imported inputs as proportion sales","Imported inputs as proprotion all inputs",
  "...of which related-party imports"), cex = .35, pt.cex = .5)

dev.off()

sum(sum3$totusinp0509)/sum(sum3$sales0509)
sum(sum3$totusinp0509)/sum(sum3$IMPANDNONIMPtotusinp0509)
sum(sum3$totusrpinp0509)/sum(sum3$IMPANDNONIMPtotusinp0509)
(sum(sum3$totusinp0509) - sum(sum3$totusrpinp0509))/sum(sum3$IMPANDNONIMPtotusinp0509)



